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I.  INTRODUCTION 


Two  different  approaches  to  the  analysis  of  a  sun  of  exponentially 
decaying  sinusoids  were  studied  In  this  report.  The  first,  and  what  is  now 
considered  the  classical  method,  was  to  minimize  the  meaui  square  difference 
between  the  data  and  the  e]q>ected  response  after  linearization  about  some 
initial  estimate  point  for  each  of  the  parameters.  This  method  was  used  most 
recently  by  Sterglopoulos^  In  the  analysis  of  decaying  Inertial  waves  excited 
In  the  fluid  contained  In  a  rotating  cylindrical  cavity  both  during  spin- up 
and  in  the  steady  state.  Details  of  this  method  are  given  In  Section  II  and 
Its  Implementation  is  described  in  Section  III.  The  second  method,  called  the 
Prony  technique,  linearizes  the  sum  of  esqponentlally  mo4Qlated  sinusoids  by 
means  of  a  recursive  definition  of  these  functions  and  thus  concentrates  the 
nonlinearity  into  a  problem  of  finding  the  complex  zeros  of  an  even  ordered 
polynomial.  Background  for  this  method  is  given  in  Section  IV  and  its 
implementation  is  given  in  Section  V. 


The  data  in  this  exercise  came  from  two  sources.  First,  several  time 
sequences  were  obtained  from  numerical  simulations  of  the  experiments  of 
Stergiopoulos  and  Aldridge.^  Results  for  this  work  are  given  in  Sections  IIIB 
and  VB  along  with  comparison  to  previous  work  by  Sedney,  Gerber,  and 
Bartos.^  A  second  source  of  data  was  from  the  gyroscope  experiments  of 
D'Amico.**  Results  of  this  work  are  given  in  Section  VIB.  All  programs  were 
left  under  the  username  AIDRIDGE  in  executable  form  so  that  any  use  of  these 
programs  should  not  require  retyping  into  a  file.  A  glossary  of  these  pro¬ 
grams  is  given  in  the  Appendix. 


II.  LINEARIZED  LEAST  SQUARES 

Here  we  give  the  details  for  the  recovery  of  frequency,  decay  rate,  amp¬ 
litude  and  phase  from  a  data  set  which  is  modeled  as  a  sum  of  exponentially 
modulated  sinusoids. 


1.  S.  Stergiopoulos t  "An  Experimental  Study  of  Inertial  Waves  in  a  Fluid 
Contained  in  a  Rotating  Cylindrical  Cavity  During  Spin-Up  From  Rest," 
Ph.D.  Thesis,  York  University,  February  1982. 

2.  S.  Stergiopoulos  and  K.D.  Aldridge,  "An  Experimental  Study  of  Complex 
Eigenfrequenoies  of  Ron-Axisymmetrio  Inertial  Waves  in  a  Rotating  Fluid 
Cylinder  During  Spin-Up  From  Rest,"  Manusoript  in  preparation. 

3.  R.  Sedney,  N.  Gerber,  and  J»M.  Bartoe,  "Oaoillations  of  a  Liquid  in  a 
Rotating  Cylinder,"  AIAA  20th  Aerospace  Sciences  Meeting,  AIAA-82-0296 , 
January  1982.  See  also  BRL  Technical  Report  ARBRL-TR-02489 ,  May  1983  (AD 
A129094). 

4.  W.P.  D'Amico,  Jr.,  W.G.  Beime,  and  T.H.  Rogers,  "Pressure  Measurements  of 
a  Rotating  Liquid  for  Impulsive  Coning  Motion,"  AIAA  20th  Aerospace 
Sciences  Meeting,  AIAA-82-0246 ,  January  1982. 
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Let  the  sequence  of  N  observations  x(n)  at  each  of  a  series  of  times 
n  A  t  be  modeled  as  a  sum  of  M  modes: 


j  «  M  -a.n 

x(n)  =•  Z  A.  sin  (u  n  +  ^  )  e  ^  +  e(n),  n  ■  1,2...N, 
j  M  1  3  3  3 

where  e(n)  is  noise  of  zero  mean  and  finite  second  moment  distribution.  Note 
that  we  do  not  assume  the  noise  is  Gaussian. 

The  parameters 

A  j .  f  and  Uj  f 

for  the  M  modes  are  unknown  and  to  be  estimated  from  the  time  sequence  x(n). 
The  calculated  response  based  on  estimates 


A , a  and  a,.  3  ^1.2. ..M 

3  3  3  3 

Is  given  by: 


j  =  M  ^  ^  ^  -an 

x(n)  =  I  A.  sin  (u  n  +  ^  )  e  ^  ,  n  *=  1,2...N. 

j  =  1  ^  J  ^ 


Our  job  is  to  estimate  the  parameters  in  a  least  squares  sense  or  otherwise. 
We  choose  to  minimize 


n  =  N  ^2 

e  =  Z  {x(n)  -  x(n) ) 
n  =  1 


where  N  is  the  number  of  observations.  If  we  minimize  with  respect  to  the 
parameters: 

A.f  w. »  and  a.,  j  «*  1,2...M 

3  j  j  3 

then  a  set  of  nonlinear  relations  would  result.  Hence  we  linearize  the 
expected  response  about  some  point  represented  by  the  superscript  o  so  that 

the  calculated  response  for  the  mode  at  time  n  becomes: 
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0^1  n)  +  (A^ 


+ 


~  ,  3x 

( u>^  -  u»  ^  )  — 


3x 

e  -  (♦i  -  ♦j)  - 


3(1) 


j 


3*, 


-  3x 

(o^  -  <*^)  — 

3  J  3a 


j 


where  the  partial  derivatives  are  evaluated  at  the  estimation  points  denoted 
by  the  superscript  o.  Summation  over  the  H  modes  gives  the  total  linearized 
response  at  time  n.  Higher  order  terms  are  ignored  and  this  will  only  lead  to 
difficulties  if  the  Initial  estimates  of  the  parameters  are  too  far  off  the 
actual  values.  Substitute  this  expression  for  x  into  the  previous  expression 
for  the  error  e  and  let 


3e  ^  3e  _  3e  _  3e  ^  ^ 

3a.  3(1).  3^^  3a. 

1  j  j  j 

This  yields  4M  equations  for  the 
obtained  from  the  Aj  equation  are 


3  *  1  >  2  .  .  .M  . 

4M  unlrnown  parameters. 


The  M  equation 


n  »  N  o 

-  E  {x(n)  -  J  (n)  -  (A 
3A.  n  »  1 
1 


3x 


3x 


3A. 


0, 


j  -  1,2...M, 


^o 

where  x  (n)  represents  the  calculated  response  at  the  superscript  o  values  of 
the  parameters.  After  collecting  terms  these  H  equations  are  written: 


n  »  N 
r 

3x 

^  i  M  -  111  > 

n  »  N 
r 

3x 

3x 

L 

n  a*  1 

3A  ^  ^ 

£• 

n  «  1 

3(i)j 

+  (♦j  -  ♦j) 


n 


E 

n  «  1 


**  3J 


3A, 


3x 

a*. 


*  -  •) 


n  «  N 

3x 

3x 

E 

n  ■  1 

3a 

3a, 

j 

j 

n  -  N  ^o 

■  E  — ^  {x(n)-x  (n))»  j  ■  1^2. ..M. 

n  •  1  3a^ 
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These  equations  can  be  written  In  matrix  form: 


T  T 

B  B  Z  »  B^R 

where  the  elements  of  the  matrices  are 


nj 


Zj  =  a  x(n)-x  (n) 


t 


n  a 


j  * 


where  p  represents  the  sequence  of  parameters  A,  to,  and  a  and  A  stands  for 
the  update  of  the  parameter.  Thus  there  are  4M  equations  In  4M  unknowns. 


Sequential  estimates  of  the  parameters 
solutions  Zj  so  that  typically 


Z^  =  A,  - 


are  obtained  from  the  update 


The  equations  are  thus  solved  recursively  until  some  convergence  criterion  is 
satisfied.  The  solution  can  be  represented  formally  as: 

Z  a  b'^^R 


It  has  been  shown  by  Hamilton^  that  if  we  let 

T  -1 

(B  B)  a  Q, 


then  the  variance  of  the  parameter  estimates  is  given  by: 


T 

R  R 

N-M  '*ii 


where  qj^j^  are  the  diagonal  elements  of  Q. 

Thus  we  obtain  error  estimates  for  each  parameter  as  well  as  the  param¬ 
eters  themselves.  The  off  diagonal  elements  give  the  covariances  of  these 
parameters.  If  the  covariances  are  near  unity  the  parameters  are  highly  cor- 

T 

related;  the  matrix  B  B  is  singular  and  there  is  no  solution.  This  is  a 


5.  W.C.  Hamilton j  Statiatioe  in  Phyeioal  Science ^  Ronald  Preea,  Neu  York, 
1964. 
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problem  of  experimental  design  in  the  sense  that  the  covariance  matrix  Is 
determined  when  the  model  Is  formulated,  which  occurs  when  the  experiment  Is 
designed.  We  shall  see  that  In  our  problem  the  exponentials  modulating  low 
frequency  sinusoids  are  themselves  not  very  orthogonal  to  each  other  so  that 
high  correlations  exist. 


III.  IMPLEMENTATION  OF  LEAST  SQUARES  METHOD 

The  algorithm  given  In  the  previous  section  was  coded  In  a  FORTRAN  pro¬ 
gram  entitled  PARTIAL. FOR  and  Is  listed  In  Appendix.  This  program  Includes 
comment  statements  so  that  It  should  be  self-explanatory.  In  order  to  supply 
confirmation  that  the  parameters  returned  by  the  program  are  Indeed  proper 
estimates  of  the  actual  amplitudes,  frequencies  and  decay  rates  some  simulated 
records  with  known  parameters  were  inverted.  The  results  of  these  simulations 
are  described  In  the  next  section. 

A.  Simulations. 


The  Immediate  application  of  this  method  was  the  Inversion  of  pressure 
data  from  numerical  experiments  by  Sedney,  Gerber  and  Barton.^  For  these 
experiments  the  frequency  and  decay  rates  were  known  approximately  for  several 
data  sets  available.  It  was  also  known  that  inversion  of  those  records  with 
low  Reynolds  numbers  was  significantly  more  difficult  than  those  at  higher 
Reynolds  numbers.  For  this  reason  a  data  set  was  simulated  with  two  modes  of 

frequencies  1.27  and  0.83  radians/second  and  decay  rates  0.45  and  0.75  sec  \ 
respectively.  In  order  to  consider  the  worst  case  as  far  u  the  Inversion  was 
concerned  both  modes  were  given  the  same  amplitude.  Results  of  Inversion  for 
this  data  are  shown  in  Table  1 .  The  parameters  of  the  Input  waveform  are 
given  above  the  horizontal  line  in  the  table;  four  cases  of  output  of  the  pro¬ 
gram,  called  SIM. FOR,  are  listed  in  the  table.  In  each  case  the  recovered  pa¬ 
rameters  are  shown  with  error  deviations  returned  by  the  Inversion  program  in 
parentheses  on  the  line  below.  The  column  labeled  "Sampling"  shows  the  number 
of  points  used  in  the  process.  The  column  titled  "Iterations"  lists  the  num¬ 
ber  of  iterations  needed  to  satisfy  the  convergence  criterion.  In  this  and  in 
most  other  cases  this  criterion  was  that  two  subsequent  outputs  differed  by 
less  than  1  part  in  1000.  Noise  was  added  to  the  last  case  and  this  will  be 
discussed  below. 

In  principle  only  9  data  points  should  be  needed  to  invert  this  data  set 
since  there  are  8  parameters.  (Phase  is  not  listed  in  the  table.)  In  pro¬ 
gressing  from  case  1  to  case  3  the  number  of  points  used  in  the  inversion  was 
decreased  from  1000  to  50.  Even  with  as  few  as  50  points  estimates  returned 
by  the  program  were  within  the  error  estimates  themselves,  which  were  still 
small.  With  fewer  points  divergence  occurs  and  this  is  due  to  round  off  error 
in  the  inversion  process.  It  might  seem  surprising  that  so  many  more  points 
than  the  minimum  were  needed  for  convergence.  Not  shown  in  this  example, 
however,  is  the  variance-covariance  matrix  of  the  inversion.  Later  it  will  be 
demonstrated  that  this  need  for  so  many  points  is  essentially  due  to  the  high 
correlation  which  exists  between  various  parameters  in  this  problem.  An  error 
in  one  produces  an  error  in  another  parameter.  Put  another  way,  decaying 
exponentials  are  not  very  orthogonal  to  one  another.  As  frequency  becomes 
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small  compared  to  decay  rate  and  ae  the  number  of  modes  in  the  Inversion 
increases,  this  situation  becomes  worse. 

The  effect  of  noise  on  the  inversion  is  shown  in  case  4.  Uniformly  dis¬ 
tributed  pseudo-random  noise  was  generated  and  ^ded  to  the  simulated  records. 
The  added  noise  had  zero  mean  and  deviation  of  10%  of  the  simulated  record. 
Even  this  modest  amount  of  noise  produces  catastrophic  results  in  the  inver¬ 
sion.  That  convergence  did  indeed  occur,  but  to  values  significantly  differ¬ 
ent  from  the  input  values,  is  at  first  somewhat  disturbing.  This  situation  is 
less  bothersome,  however,  idien  one  realizes  that  the  the  pseudo  noise  gener¬ 
ated  by  the  machine  is  indeed  correlated  to  some  extent  in  violation  of  our 
initial  assumption  of  uncorrelated  noise. 

A  more  practical  problem  for  the  inversion  of  numerical  simulations  or 
real  experimental  data  is  that  it  is  usually  not  known  how  many  modes  are 

present  in  the  record.  This  effect  was  considered  in  the  simulations  as  fol¬ 

lows.  A  record  was  simulated  with  no  added  noise  using  the  same  frequencies, 
decay  rates  and  amplitudes  as  in  Table  1 .  It  was  assumed  for  the  inversion 
that  only  one  mode  was  present  in  the  simulated  record.  Clearly  the  second 

mode  which  was  of  equal  amplitude  would  act  as  "noise"  for  the  Inversion. 

The  results  of  this  test  are  shown  in  Figures  1  and  2  which  are  graphs  of 
recovered  frequency  and  decay  rate  during  ringdown.  Each  point  is  plotted  in 
time  at  the  centre  of  a  100-point  (2-second)  interval  over  which  the  inversion 
took  place.  It  is  clear  from  the  figure  that  after  about  800  points  both  de¬ 
cay  rate  and  frequency  approach  the  initial  values  for  the  mode  recovered. 
The  interpretation  of  this  result  is  that  until  the  faster  decaying  mode  damps 
out  erroneous  values  of  both  frequency  and  decay  rate  occur.  A  practical 
approach  to  the  inversion  of  real  data  where  the  number  of  modes  is  unknown 
appears  to  be  to  track  the  inversion  process  along  the  record  in  this  manner. 
Values  of  frequency  and  decay  rate  at  the  end  of  the  record  for  the  one  mode 
recovered  can  then  be  used  in  conjunction  with  estimates  of  parameters  for  a 
second  mode.  It  was  noted  above  that  more  than  the  theoretical  minimum  number 
of  data  points  is  needed  to  invert  these  data  sets.  This  can  result  in  diver¬ 
gence  if  too  few  data  points  are  used  in  the  data  set.  It  was  found  that  no 
matter  how  slowly  the  number  of  data  points  was  decreased  from  a  number  which 
allowed  convergence  eventually  divergence  would  occur.  The  example  of  Table  1 
was  modified  in  several  ways  in  an  attempt  to  alleviate  this  problem.  The 
only  effective  change  which  allowed  fewer  data  points  without  divergence  was 
Increasing  the  frequency  of  the  modes  in  the  simulated  data.  "Hiis  clearly  had 
the  effect  of  separating  the  modes  out  to  allow  inversion  with  fewer  data 
points.  Confirmation  of  this  hypothesis  came  from  looking  at  the  variance- 
covariance  matrix  of  the  parameters  in  the  vicinity  of  divergence  due  to  small 
numbers  of  data  points.  Some  correlation  coefficients  came  within  unity  to  1 
part  in  10,000  just  before  overflow,  which  is,  of  course,  guaranteed  if  these 

T 

cofflcients  are  equal  to  1  because  the  coefficient  matrix,  B  B  ,  is  singular 
in  this  case.  What  this  means  to  the  linearized  least  squares  method  is  that 
the  signal -to-noise  ratio  must  be  increased  to  counter  the  lack  of  orthog¬ 
onality  in  cases  where  the  decay  rate  is  comparable  to  the  frequency.  This 
problem  is  inherent  in  the  type  of  function  being  used  in  the  model  which  Is^ 
of  course,  determined  by  the  experiment  being  carried  out.  The  final  solution 
to  this  problem  is  to  do  a  different  experiment  or  use  some  other  inversion 
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scheme.  It  Is  shown  later  the  the  so-called  Prony  method  has  the  potential  to 
resolve  this  problem. 

B.  Results  for  Numerical  Experiments. 

Data  from  numerical  experiments  by  Sedney,  Gerber  and  Bartos  was  Inverted 
to  find  best  estimates  of  frequency  and  decay  rate  in  the  state  of  solid  body 
rotation  for  some  axisymmetric  inertial  waves.  These  numerical  experiments 
were  carried  out  for  several  Reynolds  numbers  ranging  from  50  up  to  1000.  As 
mentioned  above,  some  difficulty  was  experienced  at  low  Reynolds  numbers  in 
resolving  modes  closely  neighboring  in  frequency.  The  data  used  in  their  ex¬ 
periments  and  also  in  this  study  irare  axial  disturbance  pressure  differences 
between  the  base  and  midplane  of  the  cylinder.  This  pressure  difference  is 
labeled  Column  3  in  the  discussion  below.  Also  available  from  this  previous 
work  was  the  pressure  difference  between  the  base  and  a  point  at  1/4  of  the 
cylinder's  height,  labeled  Column  1  below.  Finally,  a  third  field  was  con¬ 
structed  from  1/2  of  Column  3  less  Column  1  in  order  to  remove  the  response  of 
modes  with  axial  wave  number  2. 

Shown  in  Table  2  are  the  estimates  of  parameters  for  the  modes  excited  in 
the  numerical  experiments  using  the  linearized  least  squares  method.  Column  1 
in  this  table  lists  the  name  of  the  file  vihich  was  the  source  of  pressure 
data.  The  names  of  these  files  contain  two  or  three  digits  which  give  the 
perturbation  frequency  multiplied  by  100.  'Hie  first  file  is  encoded  with  the 
Reynolds  number  a  50.  The  Reynolds  number  used  in  the  numerical  e]q>eriment  is 
given  in  Column  2.  The  number  of  points  used  in  the  inversion  is  given  in 
Column  3.  Column  4  lists  the  number  of  the  pressure  field  column  21s  described 
above.  Finally,  Columns  5  through  8  give  the  estimates  of  amplitude,  frequen¬ 
cy,  phase  and  decay  rate  for  the  inversion.  Below  each  of  these  values  is  an 
error  estimate,  shown  in  parentheses,  returned  by  the  Inversion  procedure.  In 
those  entries,  with  two  rows  for  each  set  of  points,  both  modes  were  recovered 
simultaneously. 

The  results  in  this  table  divide  into  two  groups  as  a  result  of  the 
spatial  filtering  described  above.  All  entries  in  the  table  that  are  labeled 
1  and  3  In  the  fourth  column  are  modes  with  axial  wave  number  k  >  2,  while 
those  labeled  2  are  modes  of  axial  wave  number  k  «  4. 

In  general,  there  is  agreement  between  the  frequencies  and  decay  rates 
shown  in  this  table  with  those  calculated  by  SGR.  Exceptions  to  this,  of 
course,  are  those  results  which  were  not  in  common  to  SGB  and  this  report. 
This  situation  arose  when  other  modes  at  nearby  frequencies  were  found  to  be 
excited  when  forcing  took  place  at  a  particular  frequency  under  study. 

It  is  important  to  note  here,  in  comparison  of  the  results  in  this  table 

with  those  of  SGB,  what  is  meant  by  the  error  estimates  given  in  the  table. 

Implicit  In  these  estimates  is  the  integrity  of  the  model.  In  particular,  the 
noise  must  be  uncorrelated  and  all  modes  must  be  found.  If  this  is  not  true 

bias  will  be  Introduced  to  both  frequency  and  decay  rate  esUnates.  This 

phenomenon  is  Illustrated  in  the  simulations  of  Section  III-A  and  in  the 
analysis  of  gyroscope  data  in  Section  VI.  In  the  experience  of  these  calcu¬ 
lations  this  effect  is  more  inportant  in  the  decay  rate  estimates  than  in 
those  for  frequency. 
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Within  the  provisions  given  above,  frequencies  and  decay  rates  given  in 
the  table  are  in  agreement  with  those  found  by  SGB.  Exceptions  to  this  are  the 
frequencies  and  decay  rates  for  the  case  DATA125,  Reynolds  number  *  1000, 
Column  2  pressure  field.  The  reason  for  this  is  unknown  at  present.  In  order 
to  see  this  problem  in  more  detail,  the  recovery  of  the  parameters  for  this 
case  and  the  case  for  Column  3  are  described  below. 

Shown  in  Figure  3  is  the  disturbance  pressure  signal  for  the  numerical 
experiment  from  DATA125,  Column  3  pressure  field.  The  abscissa  is  the  time  in 
seconds  since  the  perturbation  stopped,  which  is  at  time  zero.  There  is  some 
indication  of  unusual  behavior  in  this  signal  near  the  point  where  t  »  8 
seconds.  Plotted  in  Figure  4  is  the  difference  between  the  pressure  signal  of 
Figure  3  and  that  calculated  at  each  point  in  time  using  the  eight  parameters 
from  the  inversion  in  Table  2.  The  large  spikes  near  t  »  8  seconds  are 
clearly  due  to  the  inability  of  the  decaying  sinusoids  to  fit  the  anomalous 
pressure  data  at  that  point.  Although  the  remaining  residual  is  only  about  2% 
of  the  maximum  pressure  signal,  it  is  important  to  note  that  this  noise  is 
rather  highly  correlated  in  its  appearance.  Furthermore,  it  is  somewhat 
nonstationary  in  that  the  statistical  properties  appear  to  be  changing  with 
time.  Both  of  these  effects  are  in  violation  of  our  assumptions  in  the 
model.  To  properly  assess  these  effects  on  the  recovered  parameters  would 
require  further  investigation. 

Data  from  Column  2  is  shown  in  Figure  5.  Figure  6  shows  the  residual 
after  only  one  mode  has  been  recovered  from  the  data.  This  plot  shows  that 
there  is  clearly  some  nonrandom  component  left  so  that  an  additional  mode  was 
sought.  This  mode  was  found  and  is  given  in  Table  1  as  having  frequency 
1 .6727  rad/sec.  “Rie  residual  shown  in  Figure  7  is  further  reduced  in  size  but 
does  display  the  same  disquieting  nonrandom,  nonstationary  noise  that  was  seen 
in  Figure  4. 

As  mentioned  above,  one  way  to  proceed  from  there  would  be  to  further 
investigate  the  effects  of  violations  of  our  assumptions.  Alternatively,  one 
could  investigate  a  noniterative  method  which  would  not  require  initial  esti¬ 
mates  of  parameters.  Then  the  uncertainty  of  the  linearized  least  squares 
method  in  knowing  whether  all  modes  has  been  found  would  be  removed.  Such  a 
method  can  be  derived  from  a  relationship  among  decaying  sinusoids  due  to 
Prony  in  1795.  We  call  this  the  Prony  Technique  and  develop  it  for  our 
problem  below. 


IV.  PRONY  TECHNIQUE 


In  the  problem  of  data  analysis  for  inertial  waves  during  ringdown  in  a 
rotating  fluid  we  seek  the  best  fit  to  data  for  modes  of  the  form: 


x(n) 


j  =  M  -a.n 

E  A.  sin  (u.n  +  ^. )  e  ^  , 
j  .  1  J  J  J 


n  ■  1,2.. .N, 


(1) 
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As  noted  earlier  in  this  report  any  direct  least  squares  procedure  for 
this  model  will  lead  to  a  nonlinear  problem  since  the  model  Is  nonlinear  In 
the  parameters  ,  4'j  >  and  . 

Our  solution  to  this  problem  was  to  linearize  about  some  arbitrary  expan¬ 
sion  point  and  solve  the  resulting  linear  least  squares  problem  Iteratively. 
This  does  Indeed  work  and  gives  parameter  estimates  with  an  error  measure.  It 
has  the  Inherent  disadvantages  of  an  Iterative  method  In  that  reasonable 
estimates  of  parameters  must  be  given  In  advance  or  divergence  might  well 
occur. 


An  alternative  procedure  to  linearization  by  expansion  as  described  above 
Is  to  linearize  by  difference  equation.  This  means  defining  our  model  recur¬ 
sively  and  in  order  to  show  this  we  represent  a  sum  of  decaying  sinusoids  as 
follows : 

* 

j  =  M  ino,  ^  -ina. 

x(n)  =  Z  A.e  ^  +  A  e  ^  n  ■  1,2...N 

j  *  1  3  J 


where  o  =  a».  +  la.  and  A  are  the  complex  eigenf requency  and  amplitude  of 
th  J  ^  ^  J 

the  j  mode,  respectively.  N  Is  the  number  of  points  and  *  means  complex 
conjugate. 

Now,  this  expression  for  x(n)  Is  linear  In  Aj  but,  of  course,  nonlinear 
in  Oj  as  mentioned  above.  As  shown  by  Chao  and  Gilbert,^  the  function  x(n) 
can  expressed  as  a  linear  combination  of  previous  values  as: 


j  =>  2M 

x(n)  »  E  S.  x(n-j)  n  =»  2M+1,...N 

j  =  1  J 


where  M  is  the  number  of  modes  used.  The  Sj  are  real  quantities  which  are 
related  to  the  through  the  following  relationship  which  Is  due  to  Prony 
and  given  in  Reference  6: 


_2M  ^  „2M-1  ^  „2M-2 
Z  -S^Z  -SjZ 


-S 


2N 


M  1 o .  la 

n  (z-e  J)  X  (z-e 


8.  B.F.  Chao  and  F.  Gilbert,  "Autoregreaaive  Eatimtion  of  Complex  Eigen- 
frequenoiea  in  Lm  Frequency  Seiamio  Spectra,"  Geophya.  J,R.  Aatro-Soc, 
£3,  841-657,  1980.  - - 
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Thus,  once  the  are  found  the  o.  can  be  found  froa  the  roots  of  the 
polynomial.  ^ 

Suppose  we  have  a  set  of  data  x(n)  and  we  wish  to  model  it  as: 

j  -  2M 

x(n)  -IS  x(n-j)  n  -  2M-t1,...N 

j  -  1  3 

where  x(n)  is  a  predicted  value  of  x(n)  based  on  the  previous  2M  values. 

Then  a  least  squares  procedure  would  be  to  minimize 


n  -  N  ^  2 

e  -  E  {x{n)  -  x(n)}  . 

n  =  2M+1 


Substitute  the  previous  expression  for  x(n)  Into  this  to  give 


n  =  N  j  -  2M  2 

£  -  E  (x{n)  -  E  S  x(n-j)}  . 

n  -  2M+1  j  -  1  ^ 


Now  to  formulate  the  least  squares  problem  we  minimize  e  with  respect  to  each 
of  the  parameters,  Sj.  Thus  setting 


g  n  =  N  j  -  2M 

“2  E  {x(n)  -  E  S  x(n-j))  x(n-k)  -0,  k  =  1,2...2M, 

k  n  -  2M+1  j  -  1  ^ 


gives  2M  equations  in  2M  unknowns.  Interchanging  the  order  of  summation  in 
the  above  relationship  gives: 


j«2M  n>»N  n>»N 

E  S  E  x(n-j)x(n-k)  ■  E  x(n}x(n-k),  k  »  1,2...2M. 

j  -  1  ^  n  -  2M+1  n  -  2M+1 


Note  that  the  terms  summed  over  n  on  both  sides  make  up  an  autocorrela¬ 
tion  of  x(n);  after  a  time  shift  the  series  is  multiplied  by  Itself  term  by 
term  and  then  added  up.  Since  this  Is  a  regression  problem  the  above  model  Is 
given  the  name  Autoregressive  or  AR  model. 
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V.  IMPLEMENTATION  PRONY  TECHNIQUE 


The  system  of  2M  linear  equations  in  2M  unknowns  is  readily  solved  to 

find  the  Sj  for  j  «  1,2..2N.  Frequency  wad  decay  rate  for  each  pair  of 

complex  roots,  z  »  exp  (lie),  can  then  be  found  from  the  earlier  polynomial 

relationship  due  to  Prony. 

A.  Simulations. 

In  order  to  test  the  Prony  method,  a  simulated  data  set  of  500  points  was 
constructed.  Tvo  sine  waves  of  frequencies  1.0  and  1.3  radians/second,  decay 
rate  0.5  1/second  and  amplitudes  0.1  and  1.0,  respectively,  were  added  togeth¬ 
er.  The  algorithm  to  recover  frequency  and  decay  rate  given  in  the  previous 
section  was  coded  as  the  program  SIMPRO.POR.  Table  3  shows  the  results  of 
this  simulation.  The  so-called  recovered  values  are  labeled  with  an  integer 
from  1  to  11.  Since  there  are  two  modes  in  the  simulated  record,  one  might 
expect  that  only  two  modes  would  need  to  be  sought.  Clearly,  %irtiat  seems  to  be 
happening,  though,  is  that  better  approximations  to  the  actual  values  of 

frequency  and  decay  rate  are  found  as  more  roots  are  sought  in  the  poly¬ 
nomial.  Finally,  with  11  roots,  the  original  values  are  reasonably  recovered. 
The  calculation  stopped  at  11  modes  sought  because  of  failure  in  the  poly¬ 
nomial  root  finder. 

It  is  important  to  note  two  points  here.  First,  no  initial  estimates 
were  required  in  carrying  out  this  calculation.  Second,  only  roots  with  posi¬ 
tive  complex  eigenfrequencies  were  retained  in  this  table;  the  program  SIMPRO 
discards  both  negative  frequencies,  decay  rates  and  any  others  which  are  out¬ 
side  the  inertial  range.  14118  last  selective  mechanism  is  somewhat  artificial 
for  the  simulated  data,  but  it  was  used  in  preparation  for  running  the  coun¬ 
terpart  of  this  simulation  program  on  real  data. 

B.  Results  for  Prony's  Method. 

The  algorithm  in  SIMPRO  was  adjusted  to  the  format  of  the  data  files  for 
the  numerical  experiments  from  BARTOS  and  called  PRONY. FOR.  A  test  of  the 
method  on  real  data  was  carried  out  on  the  BARTOS  file  OATA125,  Reynolds 
Number  1000,  438  points  and  Column  3  of  pressure.  Results  of  this  operation 
are  given  in  Table  4.  Details  of  the  input  file  DATA125  are  given  in  the  top 
of  this  table.  Not  surprisingly,  a  similar  approach  to  some  limiting  values 
of  complex  eigenf requency  is  seen  here  as  in  the  previous  simulation.  The 
last  case  is  somewhat  deviant  but  this  occurred  just  before  the  calculation 
failed  on  the  next  try,  i.e.,  9.  Hie  values  of  frequency  and  decay  rate  just 
before  this  failure  agree  within  a  few  percent  with  those  for  the  larger  mode 
in  the  last  two  rows  of  Table  2.  It  is  in^rtant  to  note  here  that  the  second 
mode  recovered  by  the  least  squares  method  and  shown  in  the  last  line  of  Table 
2  was  not  recovered  using  the  Prony  method.  This  mode  is  about  15%  of  the 
larger  mode  and  is  proportionately  larger  than  the  second  mode  used  in  the 
simulations  of  Table  3.  In  this  sense  it  might  at  first  seen  surprising  that 
the  smaller  amplitude  mode  was  not  recovered  using  the  Prony  method.  It  must 
be  remembered,  however,  that  there  was  no  added  noise  in  our  simulations  of 
Table  3.  What  is  apparently  happening  here  is  that  the  lack  of  orthogonality 
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of  these  two  modes  demands  a  larger  slgnal-to-nolse  ratio  than  is  available  in 
the  data  for  the  numerical  experiment. 


Further  work  on  the  Prony  method  seems  warranted  from  the  above  results. 
In  particular,  resolution  of  modes  should  improve  if  the  problem  were  simply 
transformed  into  the  frequency  domain,  as  in  Reference  6,  so  that  some  filter¬ 
ing  would  be  possible.  The  greatest  attraction  of  this  method  is  that  no 
initial  estimates  of  parameters  are  needed. 


VI.  FORCED  CONINS  EXPERIMEtTrS 

The  methods  developed  in  this  work  were  used  to  invert  data  from  experi¬ 
mental  measurements  made  by  D'Amico**  on  pressure  response  of  a  coning  fluid 
cylinder.  These  experiments  were  conducted  as  follows.  A  fluid  cylinder  of 
mean  radius,  a  =  3.1761  cm  and  mean  half  height,  c  =  9.9986  cm  was  set  into 
solid  body  rotation  at  a  rate  near  83.3  Hertz.  The  axis  of  symmetry  of  the 
container  which  coincided  with  the  rotation  axis  was  then  forced  to  "cone"  at 
a  fixed  angle  (0.05  degrees)  to  a  reference  irection  near  the  vertical  and  at 
various  selected  rates  near  4.0  Hertz.  After  several  seconds  the  coning  of 
the  container  was  stopped  while  the  rotation  continued.  While  the  forced 
coning  continued  and  for  several  seconds  after  it  was  stopped  pressure  on  the 
base  plate  of  the  cylinder  was  measured  and  recorded  in  analogue  form. 

The  purpose  of  the  data  analysis  carried  out  here  was  to  determine  com¬ 
plex  eigenfrequencies  from  the  freely  decaying  part  of  the  records.  Six 
analogue  records  were  digitized  and  put  into  disc  files  on  the  Vax  by  David 
Hepner  using  a  FORTRAN  program  written  by  Steve  Kushubar  for  that  purpose. 
The  Analogue  to  Digital  converter  used  in  this  operation  had  a  12-bit  resolu¬ 
tion  for  inputs  in  the  range  from  -5  volts  to  +5  volts.  The  analogue  signal 
levels  were  sufficiently  smaller  than  this  range  so  that  the  maximum  digital 
range  for  the  largest  signal  was  540  rather  than  4096.  Sampling  of  the  record 
vas  at  10,000  hertz  so  that  there  were  approximately  125  points  per  cycle  for 
signals  near  the  spin  frequency.  In  most  cases  this  was  more  than  enough  so 
that  not  every  point  was  used  in  the  data  analysis. 

Six  pairs  of  sequentially  written  files  were  provided  for  experimental 
runs  at  six  different  coning  periods.  Each  pair  of  files  consisted  of  a 
pressure  file  and  a  coning  pulse  file.  The  coning  pulse  file  consisted  of  a 
time  sequence  of  voltages  on  the  same  time  base  as  the  pressure  file.  Coning 
pulses  occurred  at  the  rate  of  one  per  coning  period  so  that  the  coning  pulse 
file  consisted  of  a  sequence  of  numbers  corresponding  to  zero  voltage  punctu¬ 
ated  by  an  integer  corresponding  to  a  voltage  pulse  on  the  original  analogue 
record  at  an  interval  of  approximately  every  2500  points. 

Each  record  consisted  of  10  seconds  of  data  so  that  there  were  100,000 
records  in  each  of  the  twelve  files  provided.  Since  sequential  files  are 
read  in  that  manner  it  is  necessary  to  read  n  records  to  get  to  the  n  1 
record.  It  would  have  taken  far  too  much  time  in  the  data  analysis  to  read 
these  files  sequentially  so  a  program  called  CONVER.FOR  was  written  to  read 
each  data  file  and  then  write  it  to  a  so-called  direct  access  file.  A  test 
program  to  look  at  the  direct  access  file  called  GETTER. FOR  showed  that  the 
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data  had  been  successfully  transferred.  This  was  necessary  because  the  direct 
access  file  could  not  be  run  by  EOT.  Haows  of  the  files  in  this  work  are 
given  in  Table  S. 

The  coning  files  were  displayed  graphically  using  the  PLOT. FOR  program  in 
order  to  find  the  location  of  the  last  coning  pulse.  This  gave  a  reference 
point  in  each  pressure  record  for  the  point  at  which  forced  aotion  had 
stopped.  Shotm  in  Column  2  of  the  above  table  is  the  approximate  point  in 
time  (1  sec  ■  10.000  points)  of  the  last  pulse.  In  two  of  the  records  there 
were  no  pulses,  which  showed  that  the  entire  pressure  record  was  teUcen  during 
ringdown  of  the  inertial  mode. 

A.  Numerical  Simulations. 


It  was  important  to  check  the  method  on  simulated  data  prior  to  running  it 
on  real  data  in  order  to  be  certain  that  the  parameters  being  recovered  were 
indeed  the  ones  in  the  original  data.  Several  data  files  were  constructed  in 
the  format  of  the  digitized  data  and  coiverted  to  direct  access  files  prior  to 
being  read  by  the  Inverting  program.  One  such  data  file  called  HEPSIM.DAT  was 
converted  to  the  direct  access  file  OIRHEP.DAT.  This  file  was  input  to  the 
Inversion  program  called  DAMICO. FOR.  Results  of  this  process  are  given  in 
Table  6.  Convergence  to  the  initial  values  is  shown  in  the  table.  The  phase 
difference  shown  is  due  to  the  fact  that  the  first  few  points  of  the  time 
sequence  were  ignored  when  the  inversion  took  place. 

B.  Results  for  Forced  Coning  Experiments. 

The  primary  response  in  the  recorded  disturbance  pressure  difference  on 
the  cylinder  baseplate  is  at  a  frequency  near  the  spin  frequency  less  the 
coning  frequency.  We  call  that  difference  the  response  frequency.  In  addi¬ 
tion  to  this  response  there  was  a  small  signal  at  the  spin  frequency  due  to 
slight  mechanical  irregularities  in  the  experimental  apparatus.  From  the 
point  of  view  of  this  data  analysis  this  signal  was  important  because  it 
allowed  us  to  determine  both  spin  and  response  during  ringdown.  both  of  which 
are  necessary  in  comparing  these  experimental  results  with  predictions  from 
theory. 

Shown  in  Table  7  are  the  recovered  parameters  from  the  six  sequences  of 
disturbance  pressure  differences  described  above.  Column  1  lists  the  names  of 
the  files  which  contained  the  pressure  data  obtained  from  the  D'Amico  experi¬ 
ments.  Column  2  lists  the  recovered  spin  rates  in  radians /second  for  five  of 
the  six  runs.  Amplitude  of  the  spin  signal  in  the  last  run  was  too  small  to 
be  successfuly  recovered.  Error  deviations  directly  from  the  inversion  pro¬ 
cess  are  given  in  parentheses  below  the  spin  values.  Recovered  frequency 
(radians/second)  and  decay  rate  (1 /seconds)  for  the  six  runs  are  shown  in 
Columns  3  and  4,  respectively.  The  columns  labeled  C^  and  C^^  are  the  frequency 

and  decay  rate  referred  to  a  fixed  (laboratory)  frame  of  reference.  The 
relationship  used  was  C^  -  1  -  response/spin  and  C^  >  decay  rate/spin. 


The  results  shown  in  Table  7  were  obtained  in  each  run  from  a  set  of  1667 
data  points  in  an  interval  during  free  decay.  The  beginning  of  the  interval 
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was  shortly  aftsr  the  last  coning  pulse  as  given  above.  Ihe  end  of  the  inter¬ 
val  was  fixed  at  the  point  which  gave  minimum  deviation  in  the  parameters  over 
a  series  of  sample  lengths.  Typically,  this  gave  records  about  2  to  2.5 
seconds  in  length. 

It  is  clear  from  the  response  and  decay  data  that  the  differences  in 
response/spin  from  one  run  to  the  next  are  at  best  only  marginally  significant 
while  the  decay/spin  differences  are  significant  enough  to  require  further 
attention.  For  if  the  response  were  explainable  by  any  linear  theory  there 
should  be  no  difference  in  these  values  with  the  different  forced  coning 
frequencies  prior  to  ringdown. 

In  order  to  assess  the  above  results  the  first  record,  PRES250,  was 
studied  in  greater  detail.  If  our  model  is  correct  then  it  should  not  matter 
what  part  of  the  Interval  during  free  decay  was  used  for  the  inversion. 
Accordingly,  each  of  the  parameters  spin,  response  and  decay  listed  in  Table  7 
were  recovered  in  subintervals  of  length  0.5  second  starting  every  0.1 
second  from  0.15  second  to  1.65  seconds  after  the  last  coning  pulse  at  0.35 
second. 

A  history  of  the  spin  during  ringdown  in  the  above  sense  is  shown  in 
Figure  8.  Each  point  is  the  recovered  spin  frequency  (radians/second)  for  a 
1667-point  Interval  centered  at  the  time  shown  on  the  abscissa.  The  vertical 
bars  are  the  standard  deviation  from  the  inversion  process.  Full  scale  on 
this  figure  is  approximately  1%  so  that  variations  in  spin  rate  during 
ringdown  appear  to  be  less  than  about  0.5%. 

Simultaneously  recovered  in  the  above  process  was  the  response  frequency 
at  each  point.  For  the  purpose  of  comparison  both  in  variability  and  devi¬ 
ation,  recovered  frequency  is  plotted  in  Figure  9  on  the  same  vertical  scale 
as  Figure  8.  Both  variation  of  the  spin  frequency  and  its  deviation  are  about 
.02%  during  ringdown.  These  small  variations  in  response  frequency  are  seen  in 
Figure  10  which  is  the  same  data  as  in  Figure  9  but  with  the  vertical  scale 
expanded  to  about  0.1%  full  scale.  The  response  frequency  is  essentially 
constant  during  ringdown. 

Had  we  not  included  the  spin  recovery  along  with  the  response  recovery  in 
these  inversions,  large  apparent  variations  in  the  observed  parameters  would 
have  occurred  during  ringdown.  This  is  best  illustrated  by  examining  re¬ 
covered  decay  rates  in  the  manner  of  subintervals  described  above  under  two 
different  conditions.  First, we  look  at  the  decay  rates  found  from  the  inver¬ 
sion  process  when  only  the  response  is  recovered.  In  this  case  the  existence 
of  the  spin  signal,  which  is  about  10%  of  the  response  at  the  beginning  of 
ringdown,  is  Ignored.  Large  variations  in  the  decay  rate  are  found  as  illus¬ 
trated  in  Figure  11.  The  deviations  in  this  figure  are  relatively  large  since 
only  500  points  were  used  in  this  example. 

For  comparison  decay  rates  found  over  the  same  interval  but  with  both 
response  and  spin  included  in  the  inversion  are  plotted  in  Figure  12  on  the 
same  vertical  scale  as  Figure  11.  The  fluctuations  evident  in  Figure  11  are 
almost  eliminated  and  an  essentially  constant  decay  rate  is  evident.  ‘Rte  mean 
value  over  the  interval  is  the  value  of  decay  rate  listed  in  Table  7  for  the 
run  PRES250. 
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The  difference  in  results  shown  in  Figures  t1  and  12,  due  to  the  spin, 
illustrates  an  important  point  regarding  errors  in  this  analysis.  With  spin 
ignored  in  the  inversion  spin  becomes  noise  as  far  as  the  inversion  is  con¬ 
cerned.  This  results  in  t%K>  effects.  First,  the  noise  is  clearly  correlated 
in  violation  of  our  assumption  of  uncorrelated,  zero  mean  noise.  Second,  all 
the  recovered  parameters  are  themselves  related  and  this  can  be  seen  in  the 
variance-covariance  matrix  of  the  parameters.  By  way  of  illustration  of  this 
correlation,  the  recovered  amplitudes  for  the  case  of  Figure  11  are  plotted 
against  recovered  decay  rates  in  Figure  13  along  with  the  errors  from  the  in¬ 
version.  Even  though  the  errors  are  large,  it  is  clear  from  the  figure  that 
amplitude  and  decay  rate  are  highly  correlated.  nils  means  that  errors  in 
calculated  amplitude  result  in  errors  in  decay  rate.  So  that  in  the  case  of 
ignoring  the  small  amplitude  spin  signal,  this  disturbance  pressure  at  this 
frequency  becomes  noise  for  the  response  at  the  free  decay  frequency.  Hence 
large  fluctuations  in  the  decay  rate  are  seen  in  Figure  11.  From  a  more 
physical  point  of  view,  it  is  clear  that  the  spin  and  response  are  separated 
by  only  about  24  rad/sec  so  that  a  "beating”  will  be  evident  in  the  combined 
effect  of  spin  plus  response.  Anomalous  values  of  decay  would  clearly  be 
found  in  such  a  case  if  a  segment  of  the  record  is  inverted  with  only  the 
response  being  modeled. 

For  comparison  with  theoretical  prediction  of  the  frequency  and  decay  rate 
data  of  Table  7,  the  last  two  columns  of  that  table  are  plotted  as  circles  in 
Figure  14.  These  data  were  scaled  by  the  values  of  and  obtained  from 

Reference  7  so  that  in  the  case  of  perfect  agreement  with  theory  all  data 
points  would  fall  at  the  point  (1,1)  shown  by  the  X  in  the  figure.  Along  with 
each  data  point  are  error  estimates  frcxn  the  inversion  process.  (Note  that 
the  vertical  and  horizontal  scales  have  been  stretched  differently  because 
frequency  is  so  much  better  determined  than  decay  rate.) 

All  the  data  for  this  (3,1,1)  mode  show  faster  decay  than  theory  predicts 
from  less  than  1%  up  to  as  much  as  65%  greater.  Frequency,  as  measured  in  the 
laboratory  frame,  is  only  slightly  but  significantly  less  than  that  predicted 
by  theory.  As  shotm  in  the  upper  corner  of  this  figure  the  Reynolds  number 
based  on  azimuthal  velocity  and  radius  of  the  cylinder  is  513,000  and  the 
half-height  to  radius  aspect  ratio  is  3.1. 

Also  shown  in  Figure  14  are  the  data  from  several  other  experiments  car¬ 
ried  out  at  York  University  by  Sterglopoulos  and  Aldridge.  These  data  include 
two  nonaxisymmetric  modes  excited  by  a  precesslng  lid  and  one  axisymmetric 
mode  excited  by  a  perturbation  in  angular  velocity.  All  of  the  latter  experi¬ 
mental  results  indicate  a  slower  decay  than  predicted.  Like  the  nonaxisym¬ 
metric  mode  studied  by  D'Amico  these  nonaxisymmetric  modes  had  slightly  lower 
frequencies  than  prediction  would  expect. 


7,  C.  V.  Kitchens j  Jr.,  N.  Gerber,  and  R.  Sedney,  "Oscillations  of  a  Liquid 
in  a  Rotating  Cylinder:  Part  I.  Solid  Body  Rotation,"  U.S.  Ballistic 
Research  Laboratory,  Aberdeen  Proving  Ground,  Maryland,  BRL  Technical 
Report  ARBRL-TR-02081 ,  June  1978  (AD  A0677S9). 
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VXX.  DXSCUSSXON 


Some  further  Interpretation  of  Figure  14  is  required.  Xt  is  noted  first 
that  the  mode  of  plotting  this  data  produces  a  differential  stretching  of  the 
errors  which  have  accumulated  in  the  experiments  and  the  data  reduction.  Those 
modes  like  the  (1,2,1)  which  have  values  of  near  zero  will  have  much  larger 

error  estimates  in  this  plot  than  those  which  are  closer  to  unity  such  as  the 
(1,1,1)  mode .  Had  we  chosen  to  ccxapare  experiment  and  theory  in  the  rotating 
frame  of  reference  this  situation  would  simply  be  reversed.  Thus  there  is  no 
best  way  to  coaqwre  but  this  differential  stretching  of  error  "bars”  must  be 
kept  in  mind  for  the  interpretation. 

Xt  is  probably  significant  that  the  best  agreement  between  theory  and  ex¬ 
periment  occurs  for  the  axisymmetric  mode.  One  is  led  to  infer  from  this  that 
the  model  used  in  the  theory  approximates  the  e]q>eriment  better  for  the 
axisymmetric  case  than  it  does  for  the  nonaxisymmebdc  case.  Xf  we  then  ask  if 
there  are  any  other  observables  in  the  nonaxisymmetric  ejqperiments  that  appear 
differently  in  the  axisymmetric  case  we  are  reminded  of  the  so-called  mean 
flow.  Xt  was  observed  by  Stergiopoulos ^  that  associated  with  the  nonaxisym¬ 
metric  mode  (1,1,1),  the  mean  flow  around  the  axis  of  rotation  became  un¬ 
stable  both  during  spln-up  from  rest  and  near  solid  body  rotation  for  very 
small  perturbations.  D'Amico  has  not  related  any  observables  to  the  mean 
flow.  Such  a  mean  flow  appears  to  become  unstable  for  axisymmetric  modes  in  a 
sphere  only  at  very  large  perturbation  amplitude.  Support  for  the  hypothesis 
that  the  mean  flow  leads  to  a  slightly  altered  spin  rate  of  the  fluid  with  a 
corresponding  shift  in  eigenfrequency  would  come  from  estimates  of  the  mean 
flow  Itself. 

Large  variations  in  the  decay  rates  are  observed  even  within  one  set  of 
experiments.  For  example, the  recovered  decay  rates  for  the  (3,1,1)  mode  vary 
by  as  much  as  65%.  The  only  experimental  variable  which  was  changed  from  one 
run  to  the  next  was  the  proximity  to  resonance  when  the  coning  was  switched 
off.  Unfortunately,  too  few  runs  have  been  processed  at  this  point  to  draw  any 
conclusions  regarding  any  possible  relationship  between  recovered  decay  rate 
and  Initial  perturbation  frequency.  More  data  will  be  processed  in  the  near 
future  to  determine  this. 
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Simulation:  2  Modes  in;  1  out 
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Figure  1.  Effect  Of  Ignoring  The  Existence  Of  A  Second  Mode  On 
Frequency  Using  Linearised  Least  Squares  Inversion. 
Points  Were  Obtained  From  100  Data  Point  Intervals 
Centered  On  The  Plotted  Location. 


Figure  2.  Same  Simulation  As  For  Fic^jre  1  Showing  The 
Effect  On  Decay  Rate. 
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Figure  3.  Axial  Disturbance  Pressure  Difference  Between  Base 

And  Center  of  Cylinder  During  Ringdown.  Abscissa  Is 
Time  Since  Perturbation  Was  Stopped  In  Seconds.  Data 
Is  From  [BARTOS]  DATA125.DAT,  c/a  =  1.000,  Reynolds 
Number  =  1000. 
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Figure  4.  Residual  After  Convergence  With  2  Modes  Modeled 
As  Listed  In  The  Last  Two  Rows  Of  Table  2. 

Same  Abscissa  As  Figure  3. 
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Figure  7.  Residual  After  2  Modes  Fitted  For  Data 

Of  Figure  5.  Parameters  For  This  Inver¬ 
sion  Are  Shown  In  Table  2,  Under  DATA125, 
1000,  400,  2. 
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Time  (seconds) 


Figure  8.  History  Of  Container  Spin  Rate  After  Coning  Was 

Stopped  For  Record  Originally  Known  as  PRES250.DAT. 
Zero  Of  Time  is  The  Start  Of  Digitization, 
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Figure  10,  Data  Of  Figure  9  Plotted  On  An  Expanded  Scale. 

Note  That  Full  Scale  In  This  Figure  Is  About  0.1%. 
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Decay  Rate:  Response  and  Spin 


D 

e 

c 

a 

y 

R 

a 

t 

e 

1 

/ 

s 

e 

c 


Time  (seconds) 


Figure  12.  Decay  Rates  During  Ringdown  For  PRES250.DAT  With 
Spin  "Noise"  Included.  Vertical  Axis  Is  The  Same 
Scale  As  In  Figure  11  For  Comparison. 
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Complex  Eigenfrequencies 
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14.  Normalized  Real  And  Imaginary  Parts  Of  Complex  Eigenfrequencies 
For  Gyroscope  Experiments  Of  D'Amico  (Circles)  And  Experiments 
Of  Stergiopoulos.  The  *  Represents  Calculated  Values  From 
Reference  7.  Legend  Is  Of  Form  [(k,n,m),  Reynolds  Number, 
Aspect  Ratio]  Where  k  Is  Axial,  n  Radial,  And  m  Azimuthal 
Wave  Number. 


TABLE  1.  EFFECTS  OF  DATA  SET  SIZE  AND  NOISE  FOR  LINEARIZED 
LEAST  SQUARES  INVERSION  OF  TWO  DECAYING  SINUSOIDS. 


An^litude 

Frequency 

Decay  Rate 

San^ling 

Iterations 

Noise 

INPUT 

1.0000 

1.27000 

0.45000 

lOOC 

1.0000 

0.83000 

0.75000 

OUTPUT 

CASE  1 

1.0000 

(0.00004) 

1.27000 

(0.00001) 

0.45000 

(0.000008) 

1000 

5 

1.0000 

(0.00007) 

0.83000 

(0.00002) 

0.75000 

(0.00004) 

CASE  2 

1.0000 

(0.00016) 

1.27000 

(0.00002) 

0.45000 

(0.00004) 

first 

300 

5 

1.0000 

(0.0001) 

0.83000 

(0.00007) 

0.75000 

(0.00017) 

CASE  3 

1.0242 

(0.19) 

1.2591 

(0.069) 

0.43915 

(0.058) 

first 

50 

12 

1.0012 

(0.019) 

0.7924 

(0.229) 

0.7283 

(0.1356) 

CASE  4 

1.5155 

(0.095) 

1.2396 

(0.032) 

0.5188 

(0.025) 

1000 

22 

10% 

0.2260 

(0.134) 

0.5117 

(0.185) 

0.3777 

(0.130) 

NOTE:  Numbers  in  parentheses  are  error  estimates. 
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TABLE  2.  INVERSION  OF  DISTURBANCE  PRESSURE  DATA  FROM  NUMERICAL 
EXPERIMENTS  OF  SGB  BY  LINEARIZED  LEAST  SQUARES. 


Bartos 

File 

Reynolds 

Number 

Points 

Col 

Amplitude 

Frequency 

Phase 

Decay 

D5083 

50 

400 

2 

-0.000514 

(0.000003) 

1.3947 

(0.0026) 

-0.582 

(0.003) 

0.4644 

(0.0035) 

DATA83B 

100 

1000 

2 

0.001587 

(0.000013) 

1.2383 

(0.0024) 

-0.410 

(0.005) 

0.2926 

(0.0030) 

DATA117 

100 

1000 

3 

0.007066 

(0.000020) 

1.1875 

(0.0011) 

0.117 

(0.002) 

0.3690 

(0.0014) 

DATA170 

100 

1000 

1 

0.003079 

(0.000129) 

0.9462 

(0.0284) 

2.847 

(0.010) 

1.2760 

(0.0229) 

DATA170 

100 

1000 

2 

-0.00143 

(0.00001) 

1.4044 

(0.0034) 

-8.600 

(0.009) 

0.2970 

(0.0029) 

DATA170 

100 

1000 

3 

-0.00310 

(0.00002) 

1.2105 

(0.0024) 

-0.151 

(0.006) 

0.2760 

(0.0021) 

NS127Z 

200 

680 

3 

0.00857 

(0.00005) 

1.3078 

(0.0006) 

2.621 

(0.004) 

0.2244 

(0.0007) 

0.00547 

(0.00005) 

0.9527 

(0.0021) 

4.791 

(0.012) 

0.4127 

(0.0035) 

NNS83Z 

200 

650 

3 

0.00382 

(0.00004) 

0.8899 

(0.0019) 

-0.723 

(0.007) 

0.3436 

(0.0025) 

-0.00217 

(0.00002) 

1.2935 

(0.0018) 

0.700 

(0.013) 

0.1801 

(0.0015) 

DATA125 

1000 

438 

1 

0.00452 

(0.00005) 

1.2909 

(0.0008) 

2.663 

(0.012) 

0.0543 

(0.0008) 

0.00153 

(0.00009) 

0.8327 

(0.0100) 

4.675 

;0.057) 

0.1438 

(0.0100) 

400 

2 

0.00337 

(0.00006) 

1.2680 

(0.0020) 

3.782 

(0.017) 

0.1181 

(0.0020) 

0.00111 

(0.00007) 

1.6727 

(0.0087) 

0.934 

(0.054) 

0.1547 

(0.0094) 

438 

3 

0.01420 

(0.00009) 

1.2819 

(0.0005) 

2.999 

(0.006) 

0.0691 

(0.0005) 

0.00244 

(0.00014) 

0.8272 

(0.0077) 

4.959 

(0.046) 

0.1310 

(0.0086) 

NOTE:  Numbers  in  parentheses  are  error  estimates. 
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TABLE  3.  INVERSION  OF  SIMULATED  DATA  USING  PHONY'S  METHOD. 

Delta  t  =  0.02000 

Number  of  points  =  500 

Actual  Values  Used  in  Simulated  Data 


Mode 

Amplitude 

Frequency 

Phase 

Decay  Rate 

1 

O.lOOOOE+01 

0.13000E+01 

O.OOOOOE+00 

0.50000E+00 

2 

O.lOOOOE+00 

O.lOOOOE+01 

O.OOOOOE+00 

0 . 50000E+00 

Recovered 

Values  of  Frequency 

and  Decay  Rate 

Modes 

Sought 

Frequency  1 

Decay  1 

Frequency  2 

Decay  2 

1 

0.12763626D+01 

0.50474519D+00 

O.OOOOOOOOD+00 

0.45665833D+00 

2 

0.12840885D+0 

0.50736495D+00 

O.OOOOOOOOD+00 

0.12690678D+02 

3 

0.12845362D4-01 

0.50813359D+00 

O.OOOOOOOOD+00 

0.55377725D+00 

4 

0.12872831D+01 

0.511 1461 5D+00 

0.38778390D+00 

0.14074909D+01 

5 

0.13004467D+01 

0.50692181D+00 

0.99295895D+00 

0.60808377D+00 

6 

0.13004877D+01 

0.50173362D+00 

0.99754570D+00 

0.52578225D+00 

7 

0.13003762D+01 

0.50081597D+00 

0.99984730D+00 

0.51301232D+00 

8 

0.13007076D+01 

0.50043101D+00 

0.10044233D+01 

0.510/4469D+00 

9 

0.13004515EH-01 

0.50022495D+00 

0.10029981D+01 

0.50632917D+00 

10 

0.13002653D+01 

0.50012187D+00 

0.10017912D+01 

0.50361231D+00 

11 

0.13000547D+01 

0. 50004 500D+00 

0. 10002 564D+01 

0.50093839D+00 
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TABLE  4.  INVERSION  OF  DATA125.DAT  USING  PRONY'S  METHOD 


Inversion  for  data  from 

nATA125: 

Delta  t 

. 

0.12500 

Shutoff  at: 

70.37500 

Frequency 

ae 

1.25000 

Number  of  points 

X 

438 

Aspect  ratio  (c/a) 

= 

1.00000 

Reynolds  number 

s 

1000.00000 

Data  column 

3 

Points  1  to  438: 

Modes  Sought 

Frequency 

Decay  Rate 

1 

0.12513086D+01 

0.46596291D+00 

2 

0.12735635D+01 

0.22293002D+00 

3 

0.12592338D+01 

0.83474899D-01 

4 

0.12597282D+01 

0.72222699D-01 

5 

0.12590036D+01 

0.67341016D-01 

6 

0.12617913D+01 

0. 6937454 5D-01 

7 

0.12620774D+01 

0.69379605D-01 

8 

0.126714630+01 

0.68961001D-01 
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TABLE  5.  NOMENCLATURE  FOR  DATA  FILES  IN  ANALYSIS  OF  GYROSCOPE  DATA. 


Coning  File 

Last  Pulse 

Hepner's  File 

Direct  Access  File 

CONE250.DAT 

3500 

PRES250.DAT 

DIRP250.DAT 

CONE255.DAT 

2000 

PRES255.DAT 

DIRP255.DAT 

CONE266.DAT 

<0 

PRES266.DAT 

DIRP266.DAT 

CONE285.DAT 

<0 

PRES285.DAT 

DIRP285.DAT 

BCON255.DAT 

5500 

BPRE255.DAT 

DIRBP255.DAT 

BCON267.DAT 

4800 

BPRE267.DAT 

DIRBP267.DAT 

TABLE  6.  DEMONSTRATION  OF  LINEARIZED  LEAST  SQUARES  INVERSION 
USING  DIRECT  ACCESS  FILE  AS  INPUT. 


Mode 

Amplitude 

Frequency 

Phase 

Decay  rate 

Starting 

estimates  are: 

1 

0.25500E+03 

0.49500E+03 

O.OOOOOE+00 

0.90000E+00 

Estimates 

after  7  iteration (s) 

0.24942E+03 

0.50000E+03 

-0.24993E+00 

0.10019E+01 

Errors 


0.86487E-01  0.13788E-02  0.34798E-03  0.13718E-02 

Actual  values  used  in  source 


250.00 


500.00 


0  1.000 
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TABLE  7.  RESULTS  FOR  INVERSION  OF  GYROSCOPE  DATA 


D '  Amico 
File 

Spin 

(Dev) 

(1/sec) 

Response 

(Dev) 

(1/sec) 

Decay 

(Dev) 

(1/sec) 

C 

r 

C. 

1 

PRESS250 

524.04 

(0.12) 

499.96 

(0.02) 

0.7797 

(0.02) 

0.04595 

0.001488 

PRESS255 

523.51 

(0.09) 

499.47 

(0.02) 

1.0953 

(0.02) 

0.04592 

0.002092 

PRESS266 

523.66 

(0.09) 

499.36 

(0.04) 

0.9480 

(0.04) 

0.04640 

0.001810 

PRESB255 

523.43 

(0.10) 

499.59 

(0.03) 

0.7430 

(0.03) 

0.04555 

0.001419 

PRESB267 

523.38 

(0.14) 

499.33 

(0.03) 

1.1807 

(0.03) 

0.04595 

0.002256 

PRESS285 

499.55 

(0.08) 

0.8841 

(0.08) 

NOTE:  Numbers  in  parentheses  are  error  estimates. 
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APPENDIX  A 


GLOSSARY  OF  PROGRAMS 


Listed  below  are  the  nanes  of  FORTRAN  programs  and  brief  descriptions  of 
their  function.  All  programs  were  left  in  executable  form  under  the  username 
AU>RIDGE,  password  XYZZY. 

In  general,  the  programs  are  self-explanatory  in  their  input  and  output. 
One  exception  to  this  is  the  use  of  CONTROL  Z  on  input.  Generally,  these  pro¬ 
grams  have  sequential  input  from  one  input  line  to  the  next.  If  one  wishes  to 
go  back  to  the  previous  input  line  instead  of  entering  a  number  on  the  current 
line,  CTRL  Z  will  usually  allow  this.  Finally,  in  the  iterative  programs, 
previously  calculated  values  of  parameters  will  be  read  into  the  input  vari¬ 
ables  by  entering  CTRL  Z  instead  of  the  requested  input.  All  inputs  in  these 
programs  are  free  format  so  that  numbers  need  only  be  separated  by  a  space. 

Numerical  Experiments 

PARTIAL. FOR  Linearized  least  squares  inversion  program  reads  input  from  BARTOS 
data  files.  Output  is  sent  to  a  file  with  extension  name  .LIS.  For  example: 
if  the  input  file  is  0ATA125.0AT  output  would  appear  in  DATA125.LIS.  Only  the 
initial  estimates  and  the  last  iteration  are  sent  to  the  .LIS  file. 

PARSEE.FOR  is  the  same  as  PARTIAL.FOR  except  that  it  displays  the  residual 
error  graphically  at  each  iteration.  This  is  useful  when  first  beginning  the 
iteration  scheme  in  order  to  display  the  data.  It  is  also  useful  after  a  mode 
has  been  recovered  to  see  what,  if  anything,  is  left. 

SIM. FOR  is  essentially  the  program  PARTIAL  set  up  for  reading  in  simulated 
data  generated  by  GEN. FOR. 

GEN. FOR  generates  1  or  more  decaying  sinusoids  for  input  to  SIM. FOR. 

Prony's  Method 

PRONY.FOR  Inverts  data  from  BARTOS  files  of  numerical  experiments  using 
Prony's  method. 

SIMPRO.FOR  used  for  testing  the  Prony  method  with  data  sets  with  known 
parameters  generated  by  SIN. FOR. 

Forced  Coning  Experiments 

DAMICO. FOR  linearized  least  squares  inversion  of  data  from  direct  access  file 
of  pressure  data  created  by  CONVER.FOR  from  sequential  data  files  of  pressure 
of  the  type  PRESS250.DAT  created  by  Hepner  for  D'Amico's  pressure  data. 

CONVER.FOR  reads  sequential  pressure  file  and  creates  a  direct  access  file 
under  the  name  DIRACC.DAT.  Operator  should  rename  this  file  appropriately. 
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GETTER. FOR  reads  direct  access  flies  created  by  CONVER.FOR  and  displays 
results  to  the  screen.  Used  for  reassurance  that  the  data  has  been  properly 
converted.  FSTDM.POR  is  the  same  program  as  DAMICO. FOR  but  without  any 
graphics. 

Miscellaneous 

KAPOINT.FOR  plots  1  or  more  sets  of  data  points  on  a  graph  and  a  legend  as 
well  as  title  and  axes.  Sample  input  for  this  marginally  documented  program 
is  CRCINM.DAT. 
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